Heatmap

Preparation

Paths and libraries setting

# Load main packages, paths and custom functions
source("../../../source/main_packages.R")
source("../../../source/functions.R")

# Florentin CONSTANCIAS's script
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_heatmap.R")

Load phyloseq objects before and after decontam

#ps <- readRDS("../../../../output/rdata/MED_phyloseq.rds")
ps.filter <- readRDS("../../../../output/1_MED/1D/1D_MED_phyloseq_decontam.rds")

Heatmap

All the MED nodes

physeq <- ps.filter
colnames(tax_table(physeq))[7] <- "Strain"

levels(physeq@sam_data$Species) <- c("AA", "CP", "CQ")
# levels(physeq@sam_data$Species)

levels(physeq@sam_data$Strain) <- c("B", "CE", "G", "L", "W-")
# levels(physeq@sam_data$Strain)

physeq %>%
  transform_sample_counts(function(x) x/sum(x) *100) %>%
  phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Strain", "Organ", "Field", "Species" ),
                          tax_aggregate = "Genus",
                          tax_add = NULL,
                          ntax  = 100) -> p
## Loading required package: ampvis2
## Warning: There are only 13 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free")+
  scale_fill_viridis_c(breaks = c(0,  0.01, 0.1, 1, 10, 100), 
                       labels = c(0,  0.01, 0.1, 1, 10, 100), 
                      trans = scales::pseudo_log_trans(sigma=0.001), # Scaling factor for the linear part of pseudo-log transformation
                       na.value = 'transparent') -> p1
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p1[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p1

Wolbachia nodes

physeq@tax_table[1, "Strain"] <- NA

physeq %>%
 transform_sample_counts(function(x) x/sum(x) *100) %>% 
 subset_taxa(Genus == "Wolbachia") %>%
 phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Strain", "Organ", "Field", "Species" ),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p
## Warning: There are only 20 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") + 
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p2
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p2[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p2

Asaia nodes

physeq %>%
 transform_sample_counts(function(x) x/sum(x) *100) %>% 
 subset_taxa(Genus == "Asaia") %>%
 phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Strain", "Organ", "Field", "Species" ),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p
## Warning: There are only 24 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") + 
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p3
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p3[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p3

Elizabethkingia nodes

physeq %>%
 transform_sample_counts(function(x) x/sum(x) *100) %>% 
 subset_taxa(Genus == "Elizabethkingia") %>%
 phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Strain", "Organ", "Field", "Species" ),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p
## Warning: There are only 5 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") + 
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p4
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p4[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p4

Legionella nodes

physeq %>%
 transform_sample_counts(function(x) x/sum(x) *100) %>% 
 subset_taxa(Genus == "Legionella") %>%
 phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Strain", "Organ", "Field", "Species" ),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p
## Warning: There are only 3 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") + 
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p5
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p5[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p5

Chryseobacterium nodes

physeq %>%
 transform_sample_counts(function(x) x/sum(x) *100) %>% 
 subset_taxa(Genus == "Chryseobacterium") %>%
 phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Strain", "Organ", "Field", "Species" ),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p
## Warning: There are only 4 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") + 
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p6
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p6[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p6

Erwinia nodes

physeq %>%
 transform_sample_counts(function(x) x/sum(x) *100) %>% 
 subset_taxa(Genus == "Erwinia") %>%
 phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Strain", "Organ", "Field", "Species" ),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p
## Warning: There are only 3 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") + 
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p7
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p7[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p7

Serratia nodes

physeq %>%
 transform_sample_counts(function(x) x/sum(x) *100) %>% 
 subset_taxa(Genus == "Serratia") %>%
 phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Strain", "Organ", "Field", "Species" ),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p
## Warning: There are only 2 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") + 
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p8
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p8[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p8

Check Blanks

# import phyloseq object before decontam
ps <- readRDS("../../../../output/1_MED/1C/1C_MED_phyloseq.rds")
colnames(tax_table(ps))[7] <- "Strain"

# 70 nodes (all)
ps %>%
  transform_sample_counts(function(x) x/sum(x) *100) %>% 
   phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Control"),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p
p + facet_grid( ~ Control, scales = "free", space = "free") +
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p9
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p9[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p9

# Pseudomonas nodes (3)
ps %>%
  transform_sample_counts(function(x) x/sum(x) *100) %>% 
  subset_taxa(Genus=="Pseudomonas") %>%
   phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Control"),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p_bis
p_bis + facet_grid( ~ Control, scales = "free", space = "free") +
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p10
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p10[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p10

# Enhydrobacter nodes (3)
ps %>%
  transform_sample_counts(function(x) x/sum(x) *100) %>% 
  subset_taxa(Genus=="Enhydrobacter") %>%
   phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Control"),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p_bis
p_bis + facet_grid( ~ Control, scales = "free", space = "free") +
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p11
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p11[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p11

# Rahnella1 nodes (3)
ps %>%
  transform_sample_counts(function(x) x/sum(x) *100) %>% 
  subset_taxa(Genus=="Rahnella1") %>%
   phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Control"),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p_bis
p_bis + facet_grid( ~ Control, scales = "free", space = "free") +
  theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p12
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p12[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p12

Save heatmaps

tiff("../../../../output/1_MED/1E/1Ed_MED_heatmap.tiff", units="in", width=28, height=16, res=300)
p1
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/1_MED/1E/1Ed_MED_heatmap_wolbachia.tiff", units="in", width=22, height=10, res=300)
p2
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/1_MED/1E/1Ed_MED_heatmap_asaia.tiff", units="in", width=22, height=10, res=300)
p3
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/1_MED/1E/1Ed_MED_heatmap_elizabethkingia.tiff", units="in", width=22, height=10, res=300)
p4
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/1_MED/1E/1Ed_MED_heatmap_legionella.tiff", units="in", width=22, height=10, res=300)
p5
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/1_MED/1E/1Ed_MED_heatmap_chryseobacterium.tiff", units="in", width=22, height=10, res=300)
p6
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/1_MED/1E/1Ed_MED_heatmap_erwinia.tiff", units="in", width=22, height=10, res=300)
p7
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/1_MED/1E/1Ed_MED_heatmap_serratia.tiff", units="in", width=22, height=10, res=300)
p8
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/1_MED/1E/1Ed_MED_heatmap_blanks_all.tiff", units="in", width=22, height=10, res=300)
p9
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/1_MED/1E/1Ed_MED_heatmap_blanks_pseudomonas.tiff", units="in", width=22, height=10, res=300)
p10
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/1_MED/1E/1Ed_MED_heatmap_blanks_enhydrobacter.tiff", units="in", width=22, height=10, res=300)
p11
dev.off()
## quartz_off_screen 
##                 2
tiff("../../../../output/1_MED/1E/1Ed_MED_heatmap_blanks_rahnella1.tiff", units="in", width=22, height=10, res=300)
p12
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/1_MED/1E/1Ed_MED_heatmap.png", units="in", width=22, height=10, res=300)
p1
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/1_MED/1E/1Ed_MED_heatmap_wolbachia.png", units="in", width=22, height=10, res=300)
p2
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/1_MED/1E/1Ed_MED_heatmap_asaia.png", units="in", width=22, height=10, res=300)
p3
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/1_MED/1E/1Ed_MED_heatmap_elizabethkingia.png", units="in", width=22, height=10, res=300)
p4
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/1_MED/1E/1Ed_MED_heatmap_legionella.png", units="in", width=22, height=10, res=300)
p5
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/1_MED/1E/1Ed_MED_heatmap_chryseobacterium.png", units="in", width=22, height=10, res=300)
p6
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/1_MED/1E/1Ed_MED_heatmap_erwinia.png", units="in", width=22, height=10, res=300)
p7
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/1_MED/1E/1Ed_MED_heatmap_serratia.png", units="in", width=22, height=10, res=300)
p8
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/1_MED/1E/1Ed_MED_heatmap_blanks_all.png", units="in", width=22, height=10, res=300)
p9
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/1_MED/1E/1Ed_MED_heatmap_blanks_pseudomonas.png", units="in", width=22, height=10, res=300)
p10
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/1_MED/1E/1Ed_MED_heatmap_blanks_enhydrobacter.png", units="in", width=22, height=10, res=300)
p11
dev.off()
## quartz_off_screen 
##                 2
png("../../../../output/1_MED/1E/1Ed_MED_heatmap_blanks_rahnella1.png", units="in", width=22, height=10, res=300)
p12
dev.off()
## quartz_off_screen 
##                 2

Session info

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS  10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ampvis2_2.6.6   forcats_0.5.0   stringr_1.4.0   dplyr_1.0.2    
##  [5] purrr_0.3.4     readr_1.4.0     tidyr_1.1.2     tibble_3.0.4   
##  [9] tidyverse_1.3.0 ggplot2_3.3.2   phyloseq_1.30.0
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-149        fs_1.5.0            lubridate_1.7.9    
##  [4] doParallel_1.0.16   RColorBrewer_1.1-2  httr_1.4.2         
##  [7] tools_3.6.3         backports_1.1.10    bslib_0.2.4        
## [10] R6_2.4.1            vegan_2.5-6         lazyeval_0.2.2     
## [13] DBI_1.1.0           BiocGenerics_0.32.0 mgcv_1.8-33        
## [16] colorspace_2.0-0    permute_0.9-5       ade4_1.7-15        
## [19] withr_2.3.0         tidyselect_1.1.0    compiler_3.6.3     
## [22] cli_2.1.0           rvest_0.3.6         Biobase_2.46.0     
## [25] network_1.16.1      xml2_1.3.2          plotly_4.9.2.1     
## [28] bookdown_0.22       sass_0.3.1          scales_1.1.1       
## [31] ggnet_0.1.0         digest_0.6.26       rmarkdown_2.7      
## [34] XVector_0.26.0      pkgconfig_2.0.3     htmltools_0.5.1.1  
## [37] dbplyr_1.4.4        htmlwidgets_1.5.2   rlang_0.4.10       
## [40] readxl_1.3.1        rstudioapi_0.11     farver_2.0.3       
## [43] jquerylib_0.1.3     generics_0.0.2      jsonlite_1.7.1     
## [46] magrittr_1.5        patchwork_1.0.1     biomformat_1.14.0  
## [49] Matrix_1.2-18       fansi_0.4.1         Rcpp_1.0.5         
## [52] munsell_0.5.0       S4Vectors_0.24.4    Rhdf5lib_1.8.0     
## [55] ape_5.4-1           lifecycle_0.2.0     stringi_1.5.3      
## [58] yaml_2.2.1          MASS_7.3-53         zlibbioc_1.32.0    
## [61] rhdf5_2.30.1        plyr_1.8.6          grid_3.6.3         
## [64] blob_1.2.1          ggrepel_0.8.2       parallel_3.6.3     
## [67] crayon_1.3.4        lattice_0.20-41     Biostrings_2.54.0  
## [70] haven_2.3.1         splines_3.6.3       multtest_2.42.0    
## [73] hms_0.5.3           ps_1.4.0            knitr_1.30         
## [76] pillar_1.4.6        igraph_1.2.6        reshape2_1.4.4     
## [79] codetools_0.2-16    stats4_3.6.3        reprex_0.3.0       
## [82] glue_1.4.2          evaluate_0.14       data.table_1.13.2  
## [85] modelr_0.1.8        vctrs_0.3.4         rmdformats_1.0.2   
## [88] foreach_1.5.1       cellranger_1.1.0    gtable_0.3.0       
## [91] assertthat_0.2.1    xfun_0.22           broom_0.7.2        
## [94] viridisLite_0.3.0   survival_3.2-7      iterators_1.0.13   
## [97] IRanges_2.20.2      cluster_2.1.0       ellipsis_0.3.1